
restore,'./local/BGPS_V1_emaf.sav',/ver
restore,'./local/BGPS_V1_pvec.sav',/ver
s = read_bgps_csv('./obsolete/bgps_distance_database.csv')

eind = where(constrain.emaf)
pvec = pvec[eind]
constrain = constrain[eind]

cnum = [4484,4357,3352]
label = ['(a)','(b)','(c)']

;; Get galactic params
mw = omni_read_conffile('./conffiles/galactic_params.conf')
dpdfs = omni_read_conffile('./conffiles/dpdf_params.conf')
conf = omni_read_conffile('./conffiles/survey_info.conf')
d = (dindgen(dpdfs.NBINS)*dpdfs.BINSIZE + dpdfs.BINSTART)/1.d3


myps,'./emaf_paper/plots/dpdf_examples_bw.eps',xsize=3.1,ysize=3.1,/bw

sun = cgSymbol('sun')

cs = 0.7

multiplot_xm,[1,3],mpcharsize=cs

FOR i=0,2 DO BEGIN
   
   j = where(emaf.cnum EQ cnum[i])
   e = emaf[j]
   k = where(pvec.cnum EQ cnum[i])
   
   ;; Build necessary arrays
   p    = pvec[k].emaf  / max(pvec[k].emaf)
   kd   = pvec[k].kdist / max(pvec[k].kdist)
   h2   = pvec[k].h2    / max(pvec[k].h2)
   dpdf = pvec[k].post  / max(pvec[k].post,jj)
   
   ;; Set titles
   xtit = i EQ 2 ? 'Heliocentric Distance  [kpc]' : ''
   ytit = i EQ 1 ? 'Relative Probability' : ''
   
   ;; Plots it all ups.
   cgPlot,d,p,charsize=cs,yr=[1d-6,1.15],/yst,thick=3,/nodata,yminor=5,$
          ytickinterval=0.5,xtit=xtit,ytit=ytit

   vline,MW.R0/1.d3*cos(e.l*!dtor)/cos(e.b*!dtor),color='blk3',thick=3,$
         linestyle=4            ; Goldenrod in the color version


   cgOplot,d,h2,thick=3,color='BLK5',linestyle=3
   cgOplot,d,kd,color='BLK5',thick=3,linestyle=2
   cgOplot,d,p,thick=3,color='blk5',linestyle=0 ; Cyan in the color version
   cgOplot,d,dpdf,thick=3,color='Black'
   
   al_legend,/top,/right,label[i],box=0,charsize=1.0*cs
   
   cgAxis,xaxis=0,xtickformat='blank_axis'
   
   ;; Mark \dml and \dbar
   cgPlots,d[jj],1.05,psym=18,symsize=0.5
   
   
   dbar = total(dpdf / total(dpdf) * d)
   cgPlots,dbar,1.075,psym=14,symsize=0.5,color='blk5' ; Magenta in the color
   
   
   
   multiplot
   
ENDFOR

myps,/done,/mp

END
